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The  addition  of  an  elevator  to  an  orbiting  space 
station  introduces  time  varying  moments  of  inertia,  a  source 
of  parametric  excitation.  This  phenomenon  of  a  parame¬ 
trically  excited  system  can  also  be  seen  when  eccentricity 
is  introduced  to  the  orbit  of  a  gravity  gradient  stabilized 
rigid  body. 

An  article  relevant  to  this  study,  "Stability  of 
Rotating  Space  Station  Containing  a  Moving  Element," 
(Salimov,  1975:41-45)  examined  the  equations  of  a  symmetric 
station  containing  two  crew  members  whose  motion  was 
described  by  either  radial  or  circular  oscillations.  The 
study  was  restricted,  however,  to  a  single  conf iguration 
vehicle  rotating  in  free  space.  Thus,  there  is  suitable 
motivation  to  consider  a  general  body  that  contains  a  moving 
element  in  an  inverse  square  gravitational  field.  That  is 
the  concern  of  the  following  study. 

I  wish  to  thank  my  thesis  advisor,  C.H.  Spenny,  of  the 
Department  of  Aeronautics  and  Astronautics,  for  his 
continuing  patience  and  assistance  throughout  the  duration 
of  this  study. 
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This  study  will  derive  the  equations  of  attitude  motion 
for  a  gravity  gradient  stabilized  space  station  whose 
moments  of  inertia  are  varying  with  time.  The  equations  are 
then  linearized,  after  which  an  analytical  solution  of  the 
pitch  equation  is  developed.  An  examination  of  the 
stability  of  motion  for  the  resulting  Hill  equation  is 
presented  and  then  compared  to  the  solution  obtained  from 
numerical  integration  of  the  nonlinear  equations.  The 
results  show  that  for  elevator  frequencies  on  the  order  of 
the  orbit  rate,  motion  can  grow  unboundedly  with  time. 
Consequently,  the  shape  of  the  classical  Lagrange  stability 
region  is  altered. 


AN  ANALYSIS  OF  SPACE  STATION  MOTION 
SUBJECT  TO  THE  PARAMETRIC  EXCITATION 
OF  PERIODIC  ELEVATOR  MOTION 

I .  Introduction 

As  the  size  of  the  space  station  and  its  platforms 
begin  to  increase,  it  will  become  necessary  to  provide  a 
means  for  travelling  from  compartment  to  compartment. 

Hence,  there  is  a  need  for  an  elevator  apparatus  that  can 
shuttle  mass  between  the  separate  modules.  Yet,  it  is  known 
that  internal  mass  movement  can  affect  the  attitude 
stability  of  a  vehicle  (Salimov,  1974:30,  1975:41). 
Consequently,  it  is  the  objective  of  this  study  to  examine 
the  motion  of  a  gravity  gradient  stabilized  space  station  in 
which  an  elevator  mechanism  operates  along  the  body  axis 
nominally  aligned  to  the  local  vertical. 

To  visualize  the  problem,  refer  to  Fig.  1.  The 
composite  center  of  mass  will  move  in  a  prescribed  circular 
orbit  with  gravity  being  the  only  external  force.  It  is 
assumed  that  the  attitude  and  orbital  motions  are  completely 
decoupled,  and  therefore  each  is  analyzed  separately.  The 
motion  of  the  elevator  relative  to  the  space  station  will  be 
prescribed  beforehand,  hence,  the  translational  motion  of 
the  elevator  and  the  space  station  relative  to  the  composite 
center  of  mass  is  known.  Thus,  with  the  orbit 


predetermined,  only  the  rotational  motion  of  the  3pace 
station  needs  to  be  analyzed. 


The  formulation  of  the  equations  of  motion  for  a  space 


station  of  arbitrary  shape  begins  with 

M  =  ‘d/dtIHl  (|) 

where  M  is  the  external  torque  produced  by  the  gravitational 
field,  and  H  is  the  angular  momentum  of  the  space  station. 
The  superscript  I  indicates  that  the  time  derivative  must  be 
taken  with  respect  to  an  inertial  frame.  The  reference 
point  about  which  to  determine  moments  and  angular  momentum 
is  the  composite  center  of  mass,  which  traces  out  the  orbit 
path.  Now,  due  to  the  prescribed  motion  of  the  elevator, 
the  "geometric  center"  of  the  space  station  moves  with 
respect  to  the  composite  center  of  mass.  The  angular 
momentum  can  be  written 


H  =  I  .'w*  (2) 

where  I  is  the  inertia  dyadic  of  the  space  station.  The 

AAA 

superscript  B  denotes  the  body  principal  basis,  bt,ba,ba, 
that  is  rigidly  attached  to  the  space  station,  while  ‘w*  is 
the  inertial  angular  velocity  of  the  space  station  with 
respect  to  an  inertial  frame.  Since  the  angular  momentum  is 
calculated  with  the  composite  center  of  mass  as  reference 


point,  the  Inertia  dyadic  I  must  be  evaluated  In  a  frame  of 
reference  not  fixed  to  the  translating  space  station  but 
attached  to  the  center  of  mass.  Thus,  I  will  be  a  function 
of  time,  I(t).  Refer  to  Eqs  (1)  and  (2),  and  permit  the 
vector  quantities  M  and  w  to  be  expressed  along  the  body 

AAA 

principal  basis,  b,,b,,ba  (See  Fig.  2).  Insert  Eq  (2)  into 
Eq  (1)  and  obtain 


M  =  »d/dt[I.,w,l  +  'w*  X  (I.'w*]  (3) 

Since  I  and  ‘w*  depend  on  time,  the  first  term  of  Eq  (3) 
becomes 


•d/dtll.'w*!  =  I  •d/dtCw" ]  +  ‘w*  d/dt[I]  (4a) 


Or,  in  matrix  form. 


•d/dtll'V]  = 


'1,0  0 

0  1,0 

0  0  I,J 

i  iw,  +  i,w, 

•  1 

=  ( I  ,Wa  +  IaWa 
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• 

1  1 
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L0  0  la. 

l  Wa 

(4b) 


(4c) 


The  second  term  of  Eq  (3)  becomes 


bi  ba  ba 


or. 


'w*  X  [i.‘s^J 


Wi  Wa  Wa 


ItWi  IaWa  IaWa 


(5a) 


!( Ia“Ia)WaWa 
(  1 1“  la  )  W1W3 
(  Ia-Ii)WiWa 


(5b) 


To  complete  Eq  (3),  consider  gravitational  torques  the  only 
input  into  M.  Using  the  inverse  square  formulation,  M  is 
given  by  (Wertz,  1978:567) 


M  =  (3u/R9)  [R  X  ( I- R  )  1  (6) 

where  R  is  the  vector  radius  from  earth  center  to  the  mass 
center  of  the  spacecraft,  and  u  =  GM  is  the  earth's 
gravitational  constant.  Note  here  that  Eq  (6)  does  not 
consider  the  gravitational  coupling  of  the  attitude  motion 
to  the  orbit,  because  terms  of  order  two  and  higher  were 
neglected  (Wertz,  1978:567).  Hence,  expand  Eq  (6)  and 
obtain  the  body  components  of  M: 

Ml  =  ( 3u/R9)  [  (  Ia-Ia)RaRa)  (7a) 

Ma  =  (  3U/R8) l ( Ii-Ia)RtRa)  (7b) 

Ma  =  (3u/Ra)t  (Ia-I.)RiRa)  (?c) 


,*•  U' 
.  '  . 

‘*JLm 

r* 
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Substituting  Eqs  (4),  (5),  and  (7)  into  Eq  (3)  yields  the 
following  equations  of  motion,  which  will  be  denoted  as 
Euler-type  equations: 

1 1W1  +  I  iWj  +  (Ii-Ia)WaWi  =  (  3u/Ra  )  (  I  a“  I  a )  RiRa  (8) 

IaWa  +  IaW,  +  (I.-Ia)WaW,  =  (  3u/R# )  (  I  a~  I »  )  RaR.  (9) 

laWa  +  IaWa  +  (Ia-Ii)WiWa  =  (  3u/R“ )  (  I  a“  1 1  )  RiRa  (10) 

When  comparing  these  equations  to  Euler's  equations  for  a 

gravity  gradient  stabilized,  rigid  body  (Kaplan,  1976:201), 
it  is  apparent  that  the  only  difference  is  the  existence  of 
the  three  Iw  terms.  These  terms  are  present  to  permit  the 
introduction  of  elevator  motion  along  arbitrary  axes. 
Moreover,  note  that  the  solution  to  Eqs  (8),  (9),  and  (10) 
will  yield  components  of  the  angular  velocity  vector  as  seen 
from  the  body  basis.  However,  the  orientation  of  the  body 
cannot  be  determined  from  these  equations.  To  obtain  the 
attitude  of  the  spacecraft,  three  angles  specifying  the 
orientation  are  introduced,  and  three  more  differential 
equations  for  these  angles  are  obtained. 

Consider  three  Euler  angles,  pitch,  yaw,  and  roll,  or 
p,  y,  and  r,  respectively,  in  a  3-1-2  series  of  Euler  angle 
rotations.  This  transformation  references  the  orientation 

AAA 

of  the  body  basis  bi/ba/ba  with  respect  to  the  orbit  basis 

A  A  A  A 

Oi,02,Oa,  where  Oi  is  the  direction  of  the  local  vertical. 


/\  /\ 

Oa  Is  tangent  to  the  orbit  path,  and  oa  Is  normal  to  the 
orbit  plane.  The  associated  direction  cosine  matrix  is 
given  by 

Aa.a(p,y,r)  =  Aa(r  )Ai(y)Aa(p) 

crcp-sysrsp  crsp+sysrcp  -cysr 
=  -cysp  cycp  sy  (11) 

srcp+sycrsp  srsp-sycrcp  cycr 

where  c  =  cosine  and  s  =  sine. 

The  expressions  for  the  rotation  angles  In  terms  of  the 
elements  of  the  direction  cosine  matrix  are 


•y. 

f-V-Vi 


p  =  -arctan(  Aai/Aaa)  (12a) 

y  =  arcsln(Aaa)  (12b) 

r  *  -arctan( Ai»/Aaa)  (12c) 

Note  that  when  y  =  90",  p  and  r  are  associated  with  the  same 
degree  of  freedom  in  the  3-1-2  sequence  of  rotations.  In 
other  words,  the  p  and  r  rotations  have  a  similar  effect  at 
y  =  90°.  (The  type  3-1-2  rotation  was  selected  because  this 
condition  of  singularity  is  unlikely  to  occur  except  for 
unstable  motion.) 

The  corresponding  Euler  rates  for  the  3-1-2  rotation 
are  given  by  (Kane  and  others,  1983:428) 


7 


1  °w,* 

pcos(y)sin(r )  +  ycos(r)j 

“Wa* 

j 

psin(y)  +  r  | 

(  °Wa* 

i 

pcos ( y) cos ( r )  +  ysin(r)] 

which  are  the  body  components  of  “w*,  the  angular  velocity 
of  the  space  station  with  respect  to  the  orbit  basis.  Thi 
angular  velocity  with  respect  to  the  inertial  basis,  ’w*  i 
needed.  Therefore, 

’w*  =  tw°  +  °v*  (14; 

And  since  the  angular  velocity  of  the  orbit  basis  with 
respect  to  the  inertial  frame,  ‘w0,  is  equal  to  Nos,  the 
eccentric  orbital  rate,  Eq  (14)  takes  the  form 


( w‘ 

An  Ata  AiaV  0 

|pcos(y)sin(r )  +  ycos(r)j 

1 Wa 

= 

Aai  Aaa  Aaa.K  0 

+|psln(y)  +  r  \ 

!  Wa 

Aai  Aaa  Aaa|[  N 

[  pcos  ( y  )cos  ( r  )  +  ysin(r)| 

w»  =  AtaN  +  pcos ( y ) s in ( r  )  +  ycos(r) 

Wa  =  AaaN  +  psin(y)  +  r 

Wa  =  AaaN  +  pcos ( y ) cos ( r )  +  ysin(r) 


The  Euler  angle  rates  can  now  be  solved  for  in  terms  of 


angular  velocity  and  Euler  angles: 

p  =  [  (W3-A33N)cos(r  )  -  (Wi-AiaN)sin(r  )  ]sec(y)  (17a) 

y  =  (Wi-A33N)cos(r  )  +  ( Wa-AasN ) s in ( r  )  (17t) 

r  =-[  (w3-Ai3n)cos(r  )  -  (Wi-Ai3N)sin(r  )  ]tan(y)  (17c) 

+  (w3-A33N) 

Since  the  torque  components  appearing  In  the  Euler-type 
equations,  Eqs  (8),  (9),  and  (10),  are  functions  of  the 
inertial  orientation  of  the  body,  they  must  be  solved 
simultaneously  with  Eqs  (17a),  (17b),  and  (17c).  This 
results  In  a  problem  consisting  of  six  first  order, 
nonlinear  differential  equations. 

Parameterization  of  the  attitude  could  also  be  carried 
out  In  terms  of  Euler  symmetric  parameters  q»,  q3,  q3,  and 
q«  instead  of  Euler  angles.  With  Euler  parameters,  there 
are  four  kinematic  equations  of  motion  that  must  be  solved 
In  conjunction  with  Eqs  (8),  (9),  and  (10).  Their 
derivation  as  well  as  the  parameterization  of  the  direction 
cosine  matrix  In  terms  of  the  quaternion  c[  is  demonstrated 
in  Appendix  A. 
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In  Eqs  (8),  (9),  and  (10),  the  Euler-type  equations  for 

a  space  station  with  elevators  operating  along  arbitrary 

•  •  • 

axes,  there  exist  the  terms  It,  la,  and  I*.  If  elevator 

• 

motion  soley  occurs  along  the  yaw  axis,  then  the  term  Ii 
vanishes.  The  significance  of  this  configuration  lies  in 
the  fact  that  Liapunov  stability  analysis  requires  the 
following  equilibrium  position  of  the  space  station:  The 
axis  of  minimum  moment  of  inertia,  or  minor  axis,  must  be 
aligned  to  the  local  vertical,  while  the  major  axis  of 
maximum  inertia  must  be  perpendicular  to  the  orbit  (Hughes, 
1986:296,298).  So,  if  the  elevator  is  to  travel  along  the 
longest  axis  of  the  space  station  inertia  ellipsoid,  then  it 

A  • 

will  operate  along  the  body  bi  direction.  Hence,  with  Ii  = 

0,  Eqs  (8),  (9),  and  (10)  become 

I iWi  +  (I3-Ia)w*w»  =  (  3u/Rs)  ( I*-Ia)RaRa  (18a) 

IaWa  +  I  iWj  4-  (I,-I,)W,W,  =  (  3  U/R9  )  (  I ,  - 1  *  )  RaR.  (18b) 

I  jWi  +  I,w»  +  ( I  a- 1 1 )  W1W2  =  (3u/R“)  (I,-I,)R,Ra  (18c) 

Now  consider  small  amplitude  motions  about  the 
equilibrium  position.  The  direction  cosine  matrix,  Eq  (11) 
becomes,  for  small  angles. 


r 


A(p,y,r)  =  -p  i  y 
r  -y  1 


Consequently,  the  radius  vector,  R,  can  be  written  in  body 
components  as 


R»  \  1  p  -r  I  R 
Ri  U  -p  1  y  |  0 
Ra  |  r  -y  1(0 


R»\-  R 


Ra  «  t-RP 


Ra  |=  \  Rr 


(21a) 

(21b) 

(21c) 


Also,  observe  that  the  orbit  bases  Oi,  oa,  o3,  rotates  at  an 
inertial  angular  velocity  N  =  Noa.  This  can  be  expressed  in 
body  coordinates 


1  p  -r  (0 

'w°  =  -p  l  y  0 

r  -y  1  In 


(22a) 


(22b) 


'  Aa  »  »  *  jn'  l 


Note  that  the  body  rates  relative  to  the  orbit  basis  can  be 
expressed  as 


•  A  9  A  »  A 

°w*  =  qbi  +  rba  +  pba. 


(23) 


When  components  are  equated  in  Eq  (14),  becomes 


Wi  = 

-rN 

• 

+  y 

(24a) 

Wa  = 

yN 

• 

+  r 

(24b) 

Wa  = 

N 

+  P 

(24c) 

Take  a  time  derivative  derivative  in  the  body  frame  of  Eqs 
(24a),  (24b)  and  (24c),  and  the  result  takes  the  form 

w,  =  y  -  (rN  +  rN)  (25a) 

w*  =  r  +  (yN  +  yN)  (2511) 

w*  =  p  +  N  (25c) 

Let  N  and  R,  for  small  eccentricity  e,  be  approximated 
(Kaplan,  1976:202) 


nil  +  2ecos(ntP)  1 

(26) 

all  -  ecos  (  nt»)  1 

(27) 

where  n 


(u/a)*'2  is  the  mean  motion,  and  t,  is  the  tii.* 


since  periapsis  passage.  (Recall  that  the  mean  anomaly  M  = 
nt„.)  Next,  differentiate  Eq  (26)  with  respect  to  time  and 
obtain 

N  =  -2n2esin(ntP)  (28) 

Finally,  substitute  Eqs  (21),  (24),  and  (25)  into  Eqs  (18a) 
(18b),  and  (18c),  and  write  the  result  as 

Ii(y-rN-Nr  )  +  (I3-Ia)  (r+yN)  ( p+N )  =  (  3u/R3)  ( Ia-Ia)  (-pR)  (rR)  (29) 
Ia(r+yN+Ny)+ia(r+yN)  +  (Ii-Ia)  (p+N) (y-rN) 

=  (3u/R3) (I, -la)  (rR) (R)  (30) 

Ia(  p+N  )+Ia(  p+N) +  (Ia-I.)  (y-rN)  (r+yN) 

=  ( 3u/Rs) ( Ia-1 1 ) ( R  )  (  -  pR  )  (31) 

Substitute  for  N  and  N  from  Eqs  (26)  and  (28),  and  delete 
the  products  of  quantities  assumed  to  be  small,  such  as  py 
and  pe .  Then,  substitute  u/R3  «  Na  into  the  above,  and  Eqs 
(29),  (30),  and  (31)  can  be  rewritten 

I>y  +  (Ia-Ia-Ii)nr  +  (Ia-I2)n2y  =  0  (32) 

lair*  +  (Ia+I.-Ia)ny  -  4(Ii-I3)n2r  +  Ia(r+yn)  =  0  (33) 
Iatp-2n2esin(nt»)  +  3n2(I2-It)p  +  I3(p  +  n)  =  0  (34) 

These  equations  are  nearly  identical  to  those  derived  by 
Kaplan  for  a  gravity  gradient  stabilized  rigid  body  (Kaplan, 


1976:204).  The  only  difference  is  the  existence  of  the  I 


terms  due  to  elevator  movement  along  the  yaw  axis.  Also, 
note  that  the  pitch  equation,  Eq  (34),  is  uncoupled  from  the 
coupled  yaw  and  roll  equations,  Eqs  (32)  and  (33). 

If  the  analysis  were  carried  out  with  elevator  motion 
along  all  three  axes,  the  character  of  the  equations  would 

remain  the  same.  This  is  because  the  IiWi  term  that  was 

•  , 

eliminated  in  Eq  (18a)  would  appear  as  I.(-rn+y)  in  Eq  (32), 


which  would  add  to  the  yaw-roll  coupling  already  present  in 
Eqs  (32)  and  (33). 


of  the 


EltCtl 


Since  the  pitch  equation  is  decoupled  from  the  roll  and 
yaw  equations,  deriving  an  analytical  solution  for  pitch 
motion  becomes  the  immediate  aim.  Hence,  insert  the 
expression  for  N  from  Eq  (26)  into  Eq  (34),  divide  by  la, 
rearrange,  and  obtain 


P  +  (  I  a/ 1  a  )  P  +  3n2[  (  Ia-I»)/Ia)p  +  (Ia/Ia)n 


=  2n*esin(ntP)  (35) 


Now,  define  the  following: 


P(t)  =  (1/2)  (1 3/ la) 

R2 ( t )  =  3nl ( Ia-I , ) /la) 


Q ( t )  =  -  ( Ia/Ia) n  +  2n2esin(nt») 


Substituting  these  expressions  into  Eq  (35)  yields  this 
linear,  second  order  equation: 


p  +  2P(t)p  +  R2  ( t ) p  =  Q ( t ) 


Introduce  a  variable  transformation,  defined  by 


Consequently,  p  can  be  written  as 

p  =  z  expl-fpdt]  (^1) 

The  necessary  time  derivatives,  p  and  *p  are 

p  =  e'/M'[z-Pzl  (^2) 

p  =  e'f'“l-2zP  +  zP1  -  zP  +  z]  (^3) 

Inserting  p,  p,  and  p  into  Eq  (35)  yields 

z  +  ( R4-P*-P  1  z  =  Qe-A-*  (W+) 

where  P,  Q,  and  R  are  given  in  Eqs  (36),  (37),  and  (38). 

The  above  is  a  linear  equation  in  which  the  coefficient  of  z 
is  a  function  of  time.  In  fact,  if  elevator  motion  is 
chosen  to  be  sinusoidal,  then  the  coefficient  of  z  is  a 
periodic  function  of  time.  Appendix  B,  accordingly,  details 
the  expansion  of  this  coefficient  as  well  as  demonstrate  the 
final  form  of  the  transformed,  homogeneous  pitch  equation  to 
be 


z"  +  [  aT+16qlcos2T+16q*cos4T  +  16qjCOs6T+16q«cos8T  ]  z  =  0 


3  ( n/w)at  ra-Ei+ ( r./2 )  ]  +  (  3/2  )  ( n/w)ar.(  r.-ra)  + 


3r  — 

( 9/8  )  (n/w)*r«2(  rj-ri-1 )  +  (13/8)r.»  -  (13/8)r.»  + 

(  27  3/128  )  r.4  +  ( 15/16 )( n/w)ar.a  (46) 

q.  =  (1/16)  [  (-3/2)  (n/w)lr.a(ra-rO  +  (  3/2 )  (n/w)ar.(ra-rt )  - 
(45/32)  (n/w)ar.a  +  (  3/2  )  ( n/w)ar.a  -  (  3/2  )  ( n/w)ar.  - 


(  93/32  )r.4  +  (  3/2  )r.a  +  r.a/2  -  r„)  (47) 

qa  =  (1/16)  [  (3/8)  (n/w)*r.a(r,-r,-l)  +  (27/32)r.4  + 

(9/16)  (n/w)ar.a  -  (1/8)  (r.a  +  (l/8)r.al  (48) 

q,  =  (1/16)  [  (-3/32)  (n/w)ar.a  +  (3/32)r.4)  (49) 

q«  *  ( 3/2048  )r,4  (50) 


and  where  n  is  the  orbit  rate,  and  w  is  the  frequency  of 
elevator  motion.  The  variables  r»  =  Ioi/Ioa  and  ra  = 

Ioa/Ioa  are  ratios  of  principal  moments  of  inertia  about  the 
center  of  mass  of  the  space  station  without  an  elevator. 
The  ratio  r0  =  MLa/4Ioa  is  the  effective  elevator  inertia 
divided  by  the  space  station  pitch  inertia. 

The  above  equation  in  z  is  a  form  of  Hill's  equation, 

z  +  (aT  +  £  gL(t) Iz  =  0  (51) 

i 

where  g'v(t)  is  a  periodic  function  (Meirovitch,  1970:280). 


This  type  of  forced  motion,  in  which  g-t(t)  acts  as  an  energy 
source,  is  an  Instance  of  parametric  excitation.  It  is 


important,  in  the  analysis  of  the  stability  of  this  motion, 
to  determine  what  values  of  aT  and  g-(t)  in  Eq  (51)  yield 
unbounded  solutions  of  the  system. 

Mathlsa.'s  Equation 

Since  the  theory  of  Hill’s  equations  has  been 
thoroughly  investigated,  the  results  can  be  immediately 
applied  to  the  parameters  of  the  present  problem.  A  special 
case  of  Hill's  equation,  Mathieu's  equation  (Hayashl, 
1964:87), 


// 

z  +  ( aT+16qcos2T ) z  =  0  (52) 

yields  a  Strutt  diagram  similar  to  the  one  shown  in  Fig.  3, 
in  which  there  is  a  partitioning  of  regions  of  stability  and 
instability  that  are  symmetric  with  respect  to  the 
horizontal  axis.  For  the  present  problem,  a  solution  is 
defined  to  be  unstable  if  it  grows  unboundedly  as  time  tends 
to  infinity.  And  on  the  other  hand,  a  solution  is  defined 
to  be  stable  if  it  remains  bounded  as  time  tends  to 
infinity.  Further,  the  variables  aT  and  q  are  functions  of 
r.,  r»,  Zi,  and  (n/w),  as  can  be  seen  from  Eqs  (46)  through 
(50).  In  other  words,  the  placement  of  points  within  the 
stable  or  unstable  zones  of  the  Strutt  diagram  depend  on  the 
vehicle  inertia  shape  (:■,  ra),  orbit  rate  to  elevator 


frequency  (n/w),  and  the  ratio  of  effective  elevator  inertia 
to  pitch  inertia  r<>. 

To  carry  out  the  stability  analysis  for  the  transformed 
pitch  equation,  one  might  use  the  fact  that  both  the  Hill 
and  Mathieu  equations  are  linear.  Therefore,  instability 
effects  in  the  Hill  equation  can  be  examined  via  Mathieu 
theory  by  inserting  the  multiple  qt  "Hill"  terms  into  the 
Mathieu  equation  one  at  a  time.  This,  of  course,  will  yield 
only  a  rough  approximation  to  the  actual  Strutt  diagram 


associated  with  Hill's  equation.  The  exact  stability  charts 


will  be  obtained  when  the  complete  nonlinear  equations  are 
numerically  integrated,  which  will  be  demonstrated  in 
Section  V. 


Mapping  ai  Efllota.  an  ths.  Strutt  diagram 

Since  the  orbiting  station  considered  here  utilizes 
only  passive,  gravity  gradient  stabilization  for  control,  it 
is  apparent  that  the  only  inertia  shapes  to  be  considered 
for  the  problem  at  hand  are  those  which  lie  in  the  classical 
Lagrange  and  Delp  regions,  as  shown  in  Fig.  4.  For  a  rigid 
body  in  a  circular  orbit,  the  Lagrange  region  is  Liapunov 
stable,  even  for  finite  motions,  whereas  the  Delp  region  is 
gyrically  stable  (Hughes,  1986:297).  Any  inertia  shape, 
i.e.,  any  combination  of  (r«,ra)  ratios  that  lies  outside  of 
these  regions  would  result  in  unstable  motion,  even  with  the 
elevator  stopped.  Moreover,  when  flexibility  effects  are 
considered  with  a  quasi-rigid  body,  the  Delp  region  loses 


-  -  r,  - 
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its  gyrlcally  stable  character  while  the  Lagrange  region 
becomes  asymptotically  three-axis  stable  (Hughes,  1986:318). 

A 

At  the  same  time,  having  an  elevator  moving  along  the  bi 
axis  o£  a  space  station  whose  inertia  shape  is  located  in 
the  Delp  region  would  correspond  to  elevator  motion  along 
the  shortest  axis  of  the  space  station  inertia  ellipsoid. 
Accordingly,  the  Delp  region  will  not  be  considered  in  this 
analysis . 

To  map  the  Lagrange  region  into  the  aT-qj  (j=l,2,3,4) 
plane,  Eqs  (46)  through  (50)  are  used.  To  carry  this  out, 
the  following  parameters  are  chosen: 

1.  r.:  0  <  r.  <  0.25 

2 .  w:  n  <_  w  <.  oo 

The  lower  limit  r0  =  0  denotes  a  null  elevator  mass.  The 
upper  limit  r<>  =  0.25  corresponds  to  the  elevator  mass  being 
equal  to  a  single  sphere  mass  of  a  dumbell  configuration 
space  station.  The  elevator  frequency  being  equal  to  the 
orbit  rate,  w  =  n,  sets  a  reasonable  lower  bound  on  elevator 
speed . 

For  each  ar-qj  diagram,  the  entire  range  of  r«  values 
are  plotted,  and  therefore  all  reasonable  elevator  masses 
are  taken  into  account.  For  each  value  of  elevator 
frequency,  w,  there  is  a  different  Strutt  diagram.  Typical 
results  are  shown  in  Fig.  5  and  Fig.  6  for  w  =  n. 

As  can  be  seen  in  Fig.  5,  the  Strutt  diagram  in  the 


ar-qt  plane  for  w  =  n,  a  polygon  shape  takes  form  beneath 
the  horizontal  axis  when  the  entire  range  of  r.  is 
considered.  In  particular/  for  a  single  r.  value,  the 
Lagrange  region  that  is  triangular  in  rt-rj  space  maps  into 
a  line  in  the  ar-qt  plane  (See  Fig.  8).  Conversely,  a  point 
in  the  two  dimensional  ar-qt  plane  will  transform  into  a 
line  segment  in  the  r»-r3  mapping,  since  several  Inertia 
shapes  give  rise  to  the  same  ar-qt  point. 

Refer  now  to  Fig.  6,  the  Strutt  diagram  in  the  aT-q* 
plane  for  w  =  n.  Observe  that  the  totality  of  feasible 
space  station-elevator  points  --  of  which  none  are  located 
in  the  unstable  region  —  lie  very  close  to  the  horizontal 
axis.  In  fact,  increasing  the  elevator  frequency,  w,  to  a 
large  value  will  condense  the  locus  of  space 
station-elevator  shapes  from  the  slender  region  shown  to  a 
spot  at  the  origin.  (This  will  be  demonstrated  later.) 

Hence,  look  now  at  the  last  two  q  variables,  q3  and  q*, 
and  notice  that  they  are  functions  of  r„,  and  not  r>  and  r2. 
The  largest  absolute  value  of  either  q3  or  q«  is  6.87  x 
10"4  when  r.  ■  0.25  and  w  =  n.  Thus,  it  is  seen  that  q3  and 
q«  are  very  small. 

On  the  whole,  take  into  account  the  facts  presented  in 
the  preceding  paragraph  as  well  as  the  actuality  of  no 
points  falling  in  the  second,  third,  or  higher  unstable 
regions  according  to  the  bounds  set  on  elevator  frequency, 
and  the  conclusion  is  made  that  only  the  first  unstable  zone 
in  the  ar-qi  plane  need  be  considered. 
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Effect  Oi.  Varying  Elevator  Frequency 

Turn  to  Fig.  7,  which  demonstrates  the  effect  of 
changing  elevator  frequency  on  the  locus  of  space 
station-elevator  points.  As  w  is  decreased  from  a  value 
equal  to  ten  times  the  orbit  rate,  lOn,  to  n,  the  first 
occurrence  of  instability  appears  at  w  =  1.9n.  Thus,  w  = 
1.9n  sets  an  approximate  ceiling  on  elevator  frequency  at 
which  unstable  motion  can  still  occur.  Or,  from  the  other 
viewpoint,  w  =  1.9n  establishes  an  approximate  lower  bound 
in  which  stable  motion  takes  place. 

Stability  Predictions 

A  determination  of  the  stability  of  the  decoupled  pitch 
motion  can  be  predicted  from  the  analytical  solution.  In 
Fig.  8,  at  r.  =  0.25,  points  1  and  7  represent  the  "fence 
posts"  enclosing  points  of  instability.  In  fact,  the  image 
of  points  1  and  7  in  the  r»-ra  plane  are  the  lines  1  and  7, 
which  also  are  boundaries  defining  the  approximate  zone  of 
instability,  as  illustrated  in  Fig.  9.  This  is  because 
points  1  and  7  enclose  the  widest  range  of  at  values,  which 
are  proportional  to  the  square  of  pitch  natural  frequency. 
Hence  any  point  lying  outside  of  the  zone  defined  by  lines  1 
and  7,  for  r.  =  0.25  and  w  =  n,  should  result  in  bounded 
pitch  motion. 

For  a  smaller  value  of  r0,  the  predicted  unstable 
region  of  Fig.  9  shrinks  in  width  as  the  lines  that  are 


labeled  line  #1  and  line  #7  for  r0 


0.25  approach  each 


other  as  r<>  Is  decreased.  The  unstable  region  will 
disappear  completely  when  r<>  =  0,  when  no  elevator  exists. 

In  Pig.  10  and  Pig.  11,  the  prediction  charts  for  an 
elevator  frequency  of  v  =  1.3n  is  shown.  When  comparing 
Fig.  9  with  Fig.  11,  -c  is  apparent  that  as  the  elevator 
frequency  is  increased,  the  location  of  the  predicted 
unstable  region  shifts.  In  fact,  as  the  elevator  frequency 
approaches  w  =  1.9n,  the  predicted  unstable  region  migrates 
toward  the  lowest  corner  of  the  Lagrange  region,  for  which 
the  only  shapes  that  result  in  unstable  motion  are  "pole” 
shaped  space  stations.  For  elevator  frequencies  greater 
than  w  =  1.9n,  linear  theory,  using  Mathieu's  equation  to 
approximate  Hill's  equation  for  the  first  unstable  region, 
predicts  no  parametrically  Induced  instability. 

Refer  to  Fig.  12  and  Fig.  13  to  see  the  effects  of  an 
elevator  frequency  w  =  (l/2)n  that's  less  that  the 
previously  set  lower  bound  of  w  =  n.  As  can  be  seen  in  Fig. 
12,  there  exist  unstable  points  in  three  unstable  regions, 
which  result  in  three  pieces  being  cut  out  of  the  Lagrange 
region  in  Fig.  13.  (It  is  important  to  note  here  that  the 
Strutt  diagram  shown  in  Fig.  12  takes  into  account  only  the 
qt  term  of  the  Hill  equation  when  plotting  the  transition 
curves.)  The  shaded  area  in  Fig.  12  represents  the  locus  of 
shapes  for  qa,  which  has  been  included  in  this  ar-qt  plot  in 
order  to  demonstrate  the  fact  that  whatever  shapes  yield  qa 
instability  will  likewise  yield  q»  instability. 

The  next  section.  Analysis  o£  tUfi.  NQhlineac  Pitch 


V.  Analysis  of  the 


Pitch  Equation 


In  the  following  cases,  an  altitude  of  150  km  in  a  low 
earth  orbit  was  chosen,  which  yields  a  circular  orbit  rate  n 
=  1.197  X  10**  sec*1  and  a  period  of  5,  250  seconds.  Several 
of  the  results  --  points  2,  4,  and  12  of  Fig.  15,  points  13, 
15,  and  20  of  Fig.  16,  points  28,  32,  and  41  of  Fig.  19,  and 
the  three  left-hand  graphs  of  Fig.  22  --  are  graphically 
displayed  only  for  a  period  of  100  orbits  in  order  to 
demonstrate  the  unique  oscillatory  shapes,  even  though  all 
points  were  run  for  a  period  of  1000  orbits.  The  remaining 
graphs  of  those  same  Figs.  15,  16,  19,  and  22  display 
results  for  a  period  of  1000  orbits  in  order  to  demonstrate 
the  unstable  character  of  motion. 

The  length  of  elevator  travel  along  the  yaw  axis  was 
chosen  to  be  90  meters.  Thus,  an  elevator  frequency  equal 
to  the  orbit  rate,  w  =  n,  corresponds  to  the  elevator 
mechanism  travelling  at  0.1143  ft/sec.  Also,  due  to  space 
limitation,  only  one  full  size  figure.  Fig.  14,  which  sets 
the  standard  for  all  the  pitch  vs.  time  graphs,  was 
included . 

Case  il:  (W  =  n) 

For  the  case  of  elevator  frequency  equal  to  the  orbit 
rate,  3  "levels"  of  r»  were  analyzed  (See  Fig.  8).  The 
results  are  shown  in  Figures  15  through  18.  Note  that  no 
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sequence  of  figures  for  r»  =  0.1  is  included  because  no 
unstable  cases  within  an  integration  time  of  1000  orbits 
could  be  located.  However,  the  sequence  of  plots  for  re  = 
0.1  seen  by  this  author  would  indicate  that  an  instability 
occurs  at  point  #24.  Also,  since  the  points  located  at  the 
level  r.  =  0.1  are  the  closest  to  the  horizontal  axis,  the 
theory  of  Mathieu  equations  indicates  that  unbounded  motion 
caused  by  e>**,  where  ji  is  the  characteristic  exponent,  would 
occur  after  a  longer  duration  of  time  since  the  magnitude  of 
p  grows  as  the  horizontal  axis  is  approached. 

In  brief,  it  can  be  easily  seen  in  Fig.  16  that  the 
linearized,  analytical  solution  gives  a  fair  approximation 
to  the  numerical  result,  which  is  shown  shaded.  The 
difference  in  shape  between  the  predicted  region  of 
instability  and  the  resulting  region  gotten  from  numerical 
integration  is  attributable  to  the  nonlinearity  of  the 
equations  of  motion,  seen  in  Eqs  (8),  (9),  and  (10). 

Further,  when  comparing  Fig.  18  to  Fig.  11,  it  is  seen  that 
the  resulting  zone  of  instability  is  "shifted  lower"  than 
the  predicted  solution. 

Case  # 2 :  (w  =  1 . 3n) 

For  elevator  frequency  equal  to  1.3  times  the  orbit 
rate,  the  results,  carried  out  for  a  single  r«  =  0.25,  are 
shown  in  Figures  19,  20,  and  21.  The  predicted  unstable 
region  (Fig.  20)  is  less  accurate  than  for  case  #1. 


Case  12:  Effect  oi. 


Motion  la  2 


One  significant  test  case  involving  roll  and  yaw  motion 
is  shown  in  Fig.  20.  The  three  left-hand  graphs  display  the 
effect  of  simple  internal  resonance  in  the  Lagrange  region, 
which  causes  bounded  exchanges  of  energy  between  pitch  and 
roll-yaw  with  a  period  of  about  30  orbits  (Breakwell  and 
Pringle,  1965:307).  For  a  value  of  r.  =  0.25  and  an 
elevator  frequency  of  w  =  n,  the  chosen  (r»,  r2)  combination 
that  produced  the  Internal  resonance  case  results  in  a  point 
"A”  (See  Fig.  18)  within  the  unstable  pitch  zone  of  the 
Strutt  diagram.  To  confirm  this  prediction,  the  elevator 
was  turned  on  at  w  =  n,  and  the  unbounded  pitch  motion 
coupled  energy  into  the  roll  and  yaw  axes  to  yield  tumbling 
along  all  three  axes. 
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This  study  has  led  to  the  following  conslusions: 

1.  For  the  problem  of  a  simple,  gravity  gradient 
stabilized,  rigid  body,  the  addition  of  elevators  introduces 
time  varying  moments  of  inertia  into  the  equations  of 
motion . 

2.  For  a  single  elevator  mechanism  whose  travel  along  the 
yaw  axis  is  described  by  a  sinusoid,  the  resulting  pitch 
motion  is  described  by  Hill's  equation,  in  which  the 
elevator  acts  as  an  energy  source  which  creates  an  instance 
of  parametric  excitation. 

3.  The  occasional  traverse  of  the  elevator  from  one  end  to 
the  other  does  not  cause  the  instability  studied  in  this 
paper.  Motion  must  be  continuous  over  many  periods. 

4.  Mathleu  stability  analysis  for  the  linearized  pitch 
equation  yields  an  approximate  lower  bound  of  elevator 
frequency  w  =  1.9n  for  stable,  bounded  pitch  motion. 

5.  The  Lagrange  region  of  stability  predicted  by  linear 
analysis  as  well  as  Liapunov's  direct  method  (Hughes, 
1986:294)  for  a  gravity  gradient  stabilized  rigid  body  may 
contain  regions  of  instability  when  an  elevator  is  operating 
along  one  of  the  principal  axes.  The  existence  of  these 
unstable  regions  require  that  the  elevator  operating  period 


be  at  or  near  the  orbit  period.  The  extent  of  the  area  of 


parametric  Instability  within  the  Lagrange  region  is 
determined  by  the  ratio  of  the  effective  moment  of  inertia 
of  the  elevator  to  the  pitch  moment  of  inertia  of  the  space 
station. 

6.  The  results  of  numerical  analysis  show  that  the 
linearized  analytical  solution  gives  a  fair  approximation  of 
the  location  of  the  1st  unstable  region. 


Recommendations 

It  is  recommended  that  an  analytical  solution  to  the 
coupled  roll-yaw  equations  be  derived,  from  which  a  basis 
can  be  formed  to  predict  the  elevator  effect  in  three 
dimensions.  One  can  expect  similar  results  in  that  further 
restrictions  be  placed  within  the  Lagrange  region.  Also,  it 
should  be  possible  to  control  disturbance  induced 
oscillations  by  making  the  elevator  position  proportional  to 
the  pitch  angle  instead  of  prescribing  the  elevator  motion 
beforehand.  Furthermore,  the  addition  of  the  elevator 
problem  might  be  considered  in  conjunction  with  the  effects 
of  flexibility  on  the  stability  of  the  space  station.  •vv"^ 


:lons  q±  Motion  using  Eular  Parana  ts  La 

Time  differentiation  of  the  quaternion  is  described  by 


d/dtlg]  =  ( 1/2 )  Zg 


(  A— 1 ) 


:512)  where 
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and  wu,  w»,  and  w.  are  the  components  of  the  angular 
velocity  V  expressed  in  the  principal  body  basis.  Since 
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Solve  for  wU/  w»,  and  w»,  and  substitute  them  into  the 
expanded  form  of  Eg  (A-l) 


And  Eg  (A-6)  can  be  rewritten  as 


q« 

q* 

q» 

q« 


(1/2)  [g2(w«-A3an)-ga(W2-A*3n)+g,(w,-A»sn)  ] 

(1/2)  ( -gJwa-Aasn) +g2(w,-AiSn) +g«(  Wi-Aaan)  1 

(1/2)  [gt(W2-A2in)-g2(Wi-Aj3n)+g«(wj-A3»n)  ] 

(1/2)  l -gt(Wi-Aun) -ga(Wa-Aaan) -ga(  Wa-Aaan)  ) 


(A-7a) 
(A-7b) 
(A -7c) 
(A-7d) 


Egs  (A-7a)  through  (A-7d)  are  the  eguivalent  set  of 
kinematic  eguations  to  be  used  instead  of  Egs  (18).  Also 
needed  in  addition  to  Egs  (A-7a)  through  (A-7d)  to  work  with 
euler  parameters  are  the  following  relations: 
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gi  =  ( 1/2 )  (gaW«-gaWv+g4Wu) 

(A-6a) 

ga  =  (1/2)  (-g»w*+gaWu+g4W„) 

• 

(A-6b) 

ga  =  ( 1/2 )  (g.w»-gawu4-g«w«) 

(A -6c) 

g<  =  (1/2)  ( -g»wu-gaw»-gaw») 

(A-6d) 
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(1/2)  ( 1  +  A»»  +  Aaa  +  Aaa+ A44 ) 


(A-8d) 


= 

which  give  euler  parameters  In  terms  of  the  elements  of  the 
direction  cosine  matrix  A.  Similarly,  A  In  terms  of  a  is 
seen  to  be 

qt2-q22-qa2+g42  2  (qiqa+q#q4)  2(qiq*-qaq4) 

A(^)  =  2(qtqa-qaq4)  -q»2+qa2-qa2+q42  2  (qaqa+q.q4)  (A-9) 

2(qiqj+qaq4)  2(qaqa-qiq4)  -qi2-qa2+qa2+q42 

Thus,  all  equations  have  been  collected  in  order  to  work 
with  euler  parameters. 


Appendix  b 


Expansion  Ql  Periodic  Coefficient 

Consider  the  dumbbell  configuration  in  Fig.  1  in  which 
the  entire  mass  of  the  space  station  minus  the  elevator 
mass,  nw,  is  concentrated  in  two  spheres  of  mass  m»,  and  the 
length  of  elevator  travel  along  the  bi  direction  is  L.  The 
moments  of  inertia  about  the  composite  center  of  mass  can  be 
described  by 


I.  =  lot  (B-la) 

la  =  1 02  +  MI  (L/2 )sin(wt )  ]*  (B-lb) 

1,  =  I oa  +  MI  (L/2  )sin(  wt )  ]a  (B-lc) 

where  M  =  [num./ (m«+m«)  ]  is  the  "reduced"  mass  of  the 
elevator,  and  I  01,  1 02,  and  I  oa  are  the  principal  moments  of 
inertia  about  the  center  of  mass  for  a  space  station  without 
the  elevator  mechanism.  A  time  derivative  of  Eq  (B-lb) 
yields 


I  =  la 

Inserting  Eq  ( B - 2 ) 

P(t) 


=  1 2  =  (wML2sin(wt )cos (wt ) ]/2 

into  Eq  (  36  )  gives  the  result 

=  fl  ( 1/2  ) wML2s in ( wt )cos ( wt ) 

2  Io2+M(  (L/2)sin(wt)  P 


(B-2) 


(B-3) 


Squaring  Eq  (B-3)  yields 

Pa(t)  =  1  I" ( 1/4  )  waMaL4slna(  wt  )cosa(wt  )~|  (B-4) 

4  [  <I0S+Mt  (L/2)sin(wt)  ]*>* 

Taking  a  time  derivative  of  Eq  (B-3)  yields 

P(t)  =  ( wML*/4  )  [  wCcos*(  wt )  -wCs  ln2(  wt )  (B -5) 

-  ( 1/ 2 )  C2s  1  na  ( wt )  cos2 ( wt )  ] 

where  C  =  [Ioa  +  ( 1/4  )MLas  ln2(  wt )  ]  ** . 

Now,  Insert  Eqs  (B-la),  (B-lb),  and  (B-lc)  Into  Eq  (37)  and 
get  the  result 

Ra(  t )  =  3naC[  Ioa+  ( 1/4  )ML*sin*(  wt )  -l0i  1  (B-6) 

Collect  all  terms,  P1,  P,  and  Ra,  from  Eqs  (B-4),  (B-5),  and 
(B-6),  respectively,  and  plug  them  into  the  transformed 
pitch  equation,  Eq  (44): 

y  +  |3naC[Io2+(l/4)MLasina(wt)-Io.]  - 

( 1/16  )CawaMaL4slna(  wt  )cosa(  wt )  -  ( 1/4  )  wMLa[Cwcosa(  wt )  - 
-  ( 1/2  )CawML*sin2(wt  )cos2(wt )  ]  j-  y 

=  Qexpl JPdt 1  (B-7) 


Cwsina(  wt ) 


Notice  the  often  occuring  term 


C  =  1 1 o*+  ( 1  /  4 ) ML2s  1  n2  ( wt )  1  "l  ^B~8) 

=  (1/Ioa  )  (1  +  x)‘‘ 

where  x  =  r«sin2(wt)  and  rB  =  ML2/4Io». 

C  can  be  expanded  using  the  binomial  series 

(l+x)“*  =  l-  x  +  x2-x*+  ....  (B-9) 

Similarly,  the  term 


C2  =  ( 1/Ioa2 )  ( 1+x ) *2 


(B-10) 


can  be  expanded  using 

(1+x)'2  =  1  -  2x  +  3xa  -  4x2  +  _  (B-ll ) 

Now,  to  decide  how  many  terms  will  approximate  each  series, 
consider  the  following:  To  determine  r.  and  Ioa,  recall  how 
the  entire  mass  of  the  rigid  space  station  without  the 
elevator  was  concentrated  in  two  spheres  of  mass  m». 
Consequently,  Ioa  =  (l/2)m«L2  and  r.  =  (l/2)(M/nu).  (This 
is  valid  since  gravity  gradient  stabilization  depends  on 
ratios  of  inertias  and  not  actual  magnitudes,  and  ratios  are 
what  are  being  expressed  here.)  Reasonable  elevator  masses 


set  a  maximum  r 


1/4,  which  corresponds  to  mi  =  m».  Also, 


performing  a  few  simple  calculations  shows  that 
approximating  (l  +  x)*‘  and  (1  +  x)'2  with  three  terms  yields  an 
inconsequential  error  while  keeping  the  number  of  terms  in 
the  approximation  to  a  minimum.  Hence,  insert  three  term 
series'  into  Eq  (B-7)  and  obtain 


z  +  J(3nVIo.)  [Ioa+(l/4)ML2sin*(wt)-Io,]  (1-x+x2) 

-  A'(l-2x+3x2)  -  B '  C '  ( 1-x  +  x2 )  +  B'D'  (l-2x  +  3x2)  j>  z  =  0  (B-12) 


where  the  homogeneous  equation  is  being  analyzed  and  the 
terms  A',  B',  C',  and  D'  are  given  by 


A'  = 


B' 

C' 

D' 


( w2M2L4/16loa2)sin*(wt  )cos2(  wt ) 
w2r,2s  i  n2  ( wt )  cos2  ( wt ) 
(W2ML2/4Ioa) 
w[cos2(  wt )  -sina(  wt )  ) 

(  WML2/2Ioa2)  s  in2  ( wt )  cos2(  wt ) 


B'C'  and  B'D'  simplifies  to 


B '  C '  =  w2r0cos(2wt) 

B'D'  =  2w2r«2s  ln2(  wt  )cos2(  wt ) 


Insert  expressions  for  B'C'  and  B'D'  into  Eq  (B-12),  carry 
out  the  various  multiplications,  and  obtain  the  result 
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z  + 


3n2r2(x) 


|3nara  +  3n2r.sln2(  wt )  -  3nari  - 
3n2r„s  lna(  wt )  (  x )  +  3nari(x)  +  3nar2(xa)  + 

3n2r<,sln2(wt )  (x2)  -  3n2rt(x2)  -  A'  +  2A'x  -  3A'xa  -  B'C' 

+  B'C’X  -  B'C'x2  +  B'D'  -  2B'D'x  *■  SB'D'x2!'  z  =  0  (B-13) 


Or, 


z  +  |3n2r2  +  3n2r.sin2(wt )  -  3n2r»  -  3n2r2r„sina(  wt )  - 

3n2r„asin4(wt )  +  3n2rir«s in2(  wt )  +•  3n2r2r.2sin4(  wt )  + 

3n2r«3sin*(  wt )  -  3n2rir„2s  in4  ( wt )  -  w2r„2sin2(  wt  )cos2(  wt ) 

+  2w2r«2s  in4(  wt  )cos2(  wt )  -  3w2r»4s  ln*(  wt  )cos2(  wt )  - 
w2r.cos  (  2wt )  +  w2r,Jcos  (  2wt )  sina(  wt )  - 
w2r.2cos(2wt)sin4(wt)  +  2w2r«2  -  4waro*s  ina(  wt ) 

+  6w2r.4sin4(  wt )  |  z  =  0  (B-14) 

Note  these  trigonometric  Identities,  and  Insert  them  Into  Eq 
( B-14 )  : 

sin2(wt )  =  1/2  -  ( 1/2 ) cos ( 2wt ) 

sin4( wt )  =  3/8  -  ( 1/2 )cos ( 2wt )  +  ( 1/8 )cos ( 4wt ) 
sln‘(wt)  =  (  9/16  )sln2(  wt )  -  (  3/16  )  cos  (  2wt ) 

+  ( 3/16 )cos ( 4wt )  -  ( l/32)cos(6wt)  +  1/32 
cos2(wt)  -  sin2(wt)  =  cos(2wt) 
cosa(wt)sina(wt)  =  1/8  -  ( 1/8  )cos  (  4wt ) 


3' 


cos2(wt)sln4(wt)  =  1/16  -  ( 1/32  )cos (  2wt ) 

-  ( 1/16 )cos ( 4wt )  +  ( 1/32 )cos ( 6wt ) 
cos2(wt)sin*(wt)  =  5/128  -  ( 1/32  )cos  (  2wt ) 

-  (1/32 )cos( 4wt )  -  ( 1/32 )cos ( 6wt )  -  ( 1/128 )cos ( 8wt ) 
sin2(wt)cos( 2wt )  =  -1/4  +  ( 1/2 )cos ( 2wt ) 

-  ( 1/4 )cos ( 4wt ) 

Eq  (B-15^  then  becomes 

z  +  |3n2rj  +  ( 1/2- ( 1/2  )  cos  (  2wt )  ]  [  3n2r.  -  3n2rar0  + 

3n2r,r.  -  4w3r«3l  -  3n*r.  +  l-3n2r.2  +  3n2r.3(ra  -  r.)  + 
6w2r.4][3/8  -  ( 1/2  )cos  (  2wt )  +  ( 1/8 ) cos ( 4wt ) ]  + 

3n2r«3[  9/32  -  (  9/32  )cos  (  2wt )  -  (  3/ 16  ) cos  (  2wt )  + 

(  3/16  )  cos  (  4wt )  -  ( 1/32  )  cos  (  6wt )  +  1/32]  -  w2r„2[l/8  - 
( 1/8  ) cos  (  4wt )  ]  +  2w2r.J[  1/16  -  ( 1/32  )cos  (  2wt )  - 
( 1/16  )cos  (  4wt )  +  (1/32  )cos  (  6wt )  -  3w*r.4[5/128  - 
( 1/32 )cos ( 2wt )  -  ( 1/32 )cos ( 4wt )  -  ( 1/32 )cos ( 6wt )  - 
( 1/128  )cos  (  8wt )  ]  -  w2r,cos(2wt)  +  w2r02[  ( 1/2  )cos  (  2wt )  - 
(1/4)  -  ( 1/4  )cos  (  4wt )  ]  -  w*r03(  (  3/8  )cos  (  2wt )  -  (1/4)  - 
( 1/4 )cos ( 4wt )  +  ( 1/16 ) cos ( 2wt )  +  ( 1/16 ) cos ( 6wt ) ] 

+  2w2r.2  |  z  =  0  (B_15) 

Group  like  terms  of  cos(2wt),  cos(4wt)/  cos(6wt),  and 
cos(8wt),  and  Eq  (B-16)  can  be  rewritten  as 

d2z/dT2  +  (aT  +  16qiCos(2T)  +  16qacos(4T)  +  16q3cos(6T) 


+  16q«cos  (  8T  )  ]  z  =  0  (B-16) 


in  which  the  Independent  variable  has  been  changed  from  t  to 
T  through  the  definition  T  =  wt,  and  the  variables  a*,  qi, 
q*,  q9,  and  q#  are  given  by 


aT  =  (n/w)a[ra-ri+(r«/2  )  ]  +  (  3/2  )  (  n/w)  2r.(  r  t-r2)  + 

(9/8)  (n/w)aro1(rJ-r,-l)  +  ( 13/8  )  r.2  -  (13/8)r.a  + 

( 27  3/128  )r.«  +  ( 15/16 )  (n/w)ar.a  (B-l?) 


q«  =  (1/16)1  (  -3/2 )  ( n/w)  ar«a(  ra-r • )  +  (  3/2 )  ( n/w)ar.( ra-rt )  - 


(  45/32)  (n/w)ar.a  +  (  3/2  )  (n/w)ar»a  -  (  3/2  )  (n/w)ar»  - 
(9  3/32) r04  +  ( 3/2 ) r.3  +  r.a/2  -  r.) 


qa  =  (1/16)1  (  3/8  )  ( n/w)  ar«a(  r*-r  »-l  >  +  (27/32)r.“  + 
( 9/16 )  ( n/w)ar«a  -  (1/8)  (r.a  +  (l/8)r.a) 


q3  =  (1/16)  [  (-3/32)  (n/w)ar.a  +  (3/32)r.*l 


q«  =  (  3/2048  )  r,* 
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Attitude  Stability  and  Natural  Frequencies  of  a  Rigid  Body 


in  a  Circular  Crbit  Plotted  Versus  r.  and  r. 
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Prediction  of  Stable  and  Unstable  Points 
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Appendix  D 


P 


Appendix  D  contains  the  Fortran  computer  program  used 
to  numerically  integrate  Eqs  (17a),  (17b),  (17c), 

( A-7a ) ,  ( A-7b) ,  ( A- 7 c ) ,  and  (A-7d). 
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£ 


8 
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program  space 

common  /ham/  t, x(7, 4), f (7,4), errest(7 ) ,n,h, 101, 102, i03,m,ll, 

+  w,e,sma,u,n4 

double  precision  t,x,f,errest,h,101,102, 103,m, ll,w,e,sma,u,n4,z, 
♦  phi, theta , ps 1 , phldot , thetadot , ps idot , phi plot, 

+  thetaplot,psiplot,a,check,conversion,n4 

dimension  z(7),  a ( 3 , 3 ) ,  phlplot(1750),  thetaplot ( 1750 ) , 

+  ps iplot ( 1750 ) 
c 
c 

c  variable  descriptions 
c 

c  t  is  time 

c  x  is  an  array  that  is  composed  of  4  copies  of  the  7  state  variables 
c  f  is  the  first  time  derivative  of  x 

c  errest  is  needed  for  the  haming  differential  equation  Integrator 
c  n  is  number  of  equations  as  well  as  the  number  of  variables 
c  h  is  the  Integration  timestep 

c  iOl,  102,  &  i 0 3  are  the  static  values  of  space  station  Inertias 


m  is  the  mass  of  the  elevator 

11  is  the  length  of  the  space  station  along  the  yaw  axis 

w  is  the  elevator  frequency 

e  is  the  eccentricity  of  the  orbit 

sma  is  the  semi-major  axis  of  the  orbit 

u  =  GM  is  the  earth's  gravitational  constant 

n4  is  the  mean  motion 

z  is  just  a  storage  place  totally  equivalent  to  x 

phi,  theta,  s  psl  are  the  euler  angles 

phldot,  thetadot,  &  psidot  are  the  euler  angle  rates 

phiplot,  thetaplot,  &  psiplot  are  plotting  storage  places 

a  is  the  direction  cosine  matrix 

check  is  a  variable  to  validate  quaternion  accuracy 


read  in  variables  needed  for  master  loops 
read  (*,*)  n,  nprint,  nstep 
read  (*,*)  t,  h 

nprint  is  the  number  of  times  data  is  printed  (outer  loop) 
nstep  is  the  number  of  steps  or  times  haming  is  called 
h  is  the  Integration  step 

total  time  of  integration  =  nprint  *  nstep  *  h 


read  in  static  inertias,  and  m,  1,  w,  e,  sma,  u,  n4 
read  (*,*)  101,  102,  103 
read  (*,*)  m,  11,  w 
read  (*,*)  e,  sma 
read  (*,*)  u 
n4  =  sqrt(u/sma**3) 


read  in  initial  conditions 
read  (*,*)  phi,  theta,  psl 
read  (*,*)  phldot,  thetadot,  psidot 


convert  euler  angles  to  radian  measure 
conversion  =  0.01745329 
phi  =  phi  *  conversion 


IeH 
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theta  =  theta  *  conversion 
psl  =  psi  *  conversion 

write  (*,*)  'phi  =  ',  phi,  'thela  =  ',  theta,  'psi  =  ',  psl 
write  (*,*)'  ' 

c 

c  compute  Initial  direction  cosine  matrix  (p5b) 

a ( 1 , 1 )  =  cos(psl)  *  cos(phl)  -  sln(theta)  *  sln(psl)  *  sln(phl) 

a ( 1 , 2 )  =  cos(psi)  *  sin(phi)  +  sin(theta)  *  sin(psi)  *  cos(phi) 

a(l,3)  =-cos(theta)  *  sin(psi) 
a(2,l)  =-cos(theta)  *  sin(phi) 
a ( 2 , 2 )  =  cos(theta)  *  cos(phi) 
a ( 2 , 3 )  =  sln(theta) 

a ( 3 , 1 )  =  sln(psi)  *  cos(phl)  +  sln(theta)  *  cos(psl)  *  sln(phi) 

a ( 3 , 2 )  =  sln(psl)  *  sln(phl)  -  sln(theta)  *  cos(psl)  *  cos(phl) 

a(3,3)  =  cos(theta)  *  cos(psl) 
c 

c  compute  the  initial  angular  velocity  (  z(5),  z(6),  z ( 7 )  are  wl, 
c  w2,  and  w3,  respectively  ) 

2(5)  =  a(l,3)*n4  +  phidot*cos(theta)*sin(psi )  +  thetadot*cos(psi ) 
z ( 6 )  =  a(2,3)*n4  *■  phldot*sln(theta)  +  pstdot 

2(7)  =  a(3,3)*n4  +  phidot*cos(theta)*cos(psi )  +  thetadot*sin(psi ) 
c 

c  compute  initial  quaternion  (p5c) 

2(4)  =  0.5  *  sqrt(  1  +  a ( 1, 1 )  +  a(2,2)  +  a(3,3)  ) 

z(l)  =  ( 0 . 25/z ( 4 ) )  *  (  a(2,3)  -  a(3,2)  ) 

z(2)  =  (0.25/z(4))  *  (  a(3,l)  -  a(l,3)  ) 

2(3)  =  (0.25/z(4) )  *  (  a ( 1, 2 )  -  a(2,l)  ) 

c 

c  check  sum  of  quaternion  components  squared  being  equal  to  one 
check  =  z(l)**2  +  z(2)**2  *■  z(3)**2  t  z(4)**2 
c 

c  print  tabular  headings  and  initial  values 
write  (*,*)  '  ' 

write  (*,*)  '  ' 

write  (*,11)  t,  (z(j),  j =5 , 7 ) ,  check,  phi,  theta,  psi,  phidot, 

+  thetadot,  psldot 

11  format  ( f 9 . 1 ,  lx,  10(f9.6,lx)) 
c 

c  rename  variables  for  hamlng-sultable  iteration 
x( 1, 1 )  =  z( 1) 
x( 2, 1 )  =  z ( 2 ) 
x( 3, 1)  =  z ( 3 ) 
x( 4, 1 )  =  z( 4) 
x( 5, 1 )  =  2(5) 
x( 6, 1 )  =  z(6) 
x( 7, 1 )  =  z ( 7 ) 
c 

c  Initialize  hamlng 
nxt  =  0 

call  haming(nxt) 

If (nxt  .ne.  0)  go  to  50 
1  format(2x, 'hamlng  did  not  initialize') 
write  (*,1) 
stop 
c 
c 

50  continue 
c 

c  integrate  ordinary  differential  equations  ...  two  nested  loops 


do  200  lpc  =  1,  nprint 
c 

do  100  istp  =  1,  nstep 
c 

call  haming(nxt) 

c  each  call  to  haming  advances  one  step 

c 

100  continue 

c  after  nstep  integration  steps,  convert  to  print  variables  and 

c  print  the  current  values 

c 

z ( 1 )  =  x(l,nxt) 

2(2)  =  x(2,nxt) 

2(3)  =  x(3,nxt) 
z ( 4 )  =  x(4,nxt) 
z(5)  =  x(5,nxt) 

2(6)  =  x(6,nxt) 

2(7)  =  x(7,nxt) 
c 

c  compute  dir  cos  matrix  a,  euler  angles,  &  check  quaternion  (p5) 
c  also,  store  data  for  plotter 

a ( 1, 3 )  =  2  *  (  2( 1 ) *2( 3 )  -  2( 2 ) *2{ 4 )  ) 
a ( 2, 1 )  =  2  *  (  2 ( 1 ) *2 ( 2 )  -  z( 3 ) *2( 4 )  ) 
a ( 2 , 2 )  =  -z(l)**2  +  2(2 )  **2  -  2(3 ) **2  +  2(4  ) **2 
a(2,3)  =  2  *  (  z( 2 ) *z( 3 )  f  2(l)*z(4)  ) 
a ( 3, 3 )  =  -2 ( 1 )  **2  -  2(2) **2  +  z(3)**2  +  2 ( 4 ) **2 
c 

phi  =  atan(  -a ( 2 , 1 ) /a ( 2 , 2 )  ) 
theta  =  asin(  a ( 2 , 3 )  ) 
psi  =  atan(  -a ( 1, 3 ) /a ( 3, 3 )  ) 
c 

phlplot(lpr)  =  phl/conversion 
thetaplot( Vpt)  »  theta/conversion 
psiplot(lpr)  =  psi/converslon 
c 

phldot  =  (  (z( 7 )-a( 3, 3) *n4 )*cos(psi )  -  (2( 5)-a(l, 3 ) *n4 ) * 

+  sln(psl)  )  /  cos(theta) 

thetadot=(  ( z( 5)-a ( 1, 3 ) *n4 ) *cos(psi )  +  (z(7)-a(3,3)*n4)* 

+  sln(psi)  ) 

psldot  =  (z(6)-a(2,3)*n4)  - 

t  (  (z(7)-a(3,3)*n4) *cos (psi )  -  (z(5)-a(l,3)*n4)* 

+  sln(psl)  )  *  tan(theta) 

c 

check  =  z(l)**2  +  z(2)**2  +  z(3)**2  +  z(4)**2 
c 

write  (*,11)  t,  (z(J),  J=5, 7 ) ,  check,  phi,  theta,  psi,  phldot, 

♦  thetadot,  psidot 

c 
c 
c 

200  continue 

c 

c 

c  output  the  euler  angles  so  that  an  edited  output  file  can  be  sent 
c  to  the  plotter 

c 

do  60  i  =  1,  nprint 

t  =  l  *  (h*nstep) 

write  (*,85)  t,  ph 1  plot ( i ) 

60  continue 


do  70  1  =  1,  nprlnt 

t  =  i  *  (h*nstep) 

write  (*,85)  t,  thetaplot(i) 

70  continue 
c 

do  80  i  =  1,  nprlnt 

t  =  1  *  (h*nstep) 

write  (*,85)  t,  psiplot(i) 

80  continue 
c 

85  format  (£9.1,  3x,  f20.13) 
stop 
end 
c 
c 
c 
c 
c 

subroutine  hamlng(nxt) 

common  /ham/  x,y(7, 4) , f (7,4) ,errest(7),n,h, 101, 102, 103, m, 11, 

+  w,e,sma,u,n4 

double  precision  x,y, f ,errest,h,xo, tol,hh,  101,  102, 103,m, 11, 

+  w,e,sma,u,n4 

tol  =  1.0d-08 
if(nxt)  190,10,200 
10  xo  =  x 

hh  =  h/2 .d+00 
call  rhs(l) 
do  40  1  =  2,4 
x  =  x  +  hh 
do  20  1  ■  l,n 

20  y ( i , 1 )  *  y( i ,1-1)  +  hh*f (1,1-1) 
call  rhs(l) 
x  =  x  +  hh 
do  30  i  =  l,n 

30  y( 1,1)  =  y(l, 1-1)  +  h*f ( 1,1) 

40  call  rhs(l) 
jsw  =  -10 
50  lsw  =  1 

do  120  1  =  l,n 

hh  =  y( 1,1)  +  h*(  9 -d  +  00*f ( 1,1)  +  19  .d  +  00*f ( 1,2)  -  5.d  +  00*f ( 1 , 3 ) 
1  4  £(1,4)  )  /  24.d400 

1  £ (  dabs(  hh  -  y ( 1 , 2 ) )  .It.  tol  )  go  to  70 
lsw  =  0 

70  y(l, 2)  =  hh 

hh  =  y ( 1,1)  4  h* (  £(1,1)  4  4 .d400*£ ( 1,2)  4  £ ( 1 , 3 ) )/3.d400 
1  £ (  dabs(  hh-y( 1,3))  .It.  tol  )  go  to  90 
lsw  =  0 

90  yd, 3)  =  hh 

hh  =  y ( 1,1)  4  h* (  3.d400*£(i,l)  4  9.d400*f(i,2)  4  9 .d400*£ ( 1,3) 

1  4  3.d400*f(i,4)  )  /  8.d400 

If (  dabs(hh-y( 1,4))  .It.  tol  )  go  to  110 

lsw  =  0 

110  yd, 4)  =  hh 
120  continue 
x  =  xo 

do  130  1  =  2,4 
x  =  x  4  h 
130  call  rhs(l) 


if(isw)  140,140,150 
jsw  =  Jsw  4  1 
lf(jsw)  50,280,280 

X  =  xo 

lsw  =  1 
jsw  =  1 

do  160  1  =  l,n 
errest(l)  =  0.0 
nxt  =  1 
go  to  280 
jsw  =  2 

nxt  =  labs(nxt) 
x  =  x  +  h 

npl  =  mod (nxt, 4)  +  1 
go  to  (210, 230), lsw 
go  to  (270, 270, 270, 220), nxt 
isw  =  2 

nm2  =  mod (npl, 4)  +  1 
nml  =  mod(nm2,4)  +  1 
npo  =  mod (nml, 4)  +  1 
do  240  1  =  l,n 

f ( i,nm2)  =  y(i,npl)  +  4.d+00*h*(  2.d+00‘f ( l,npo)  -  f(i,nml) 

1  ♦  2.d+00*£(l,nm2)  )  /  3.d+00 

y(l,npl)  =  £ ( 1 ,nm2 )  -  0.925619835d+00*errest( 1 ) 
call  rhs(npl) 
do  250  i  =  l,n 

y(l,npl)  =  (  9.d+00*y(i,npo)  -  y(l,nm2)  +  3.d+00*h*(  £(l,npl) 
1  +  2 .d+00*£( l,npo)  -  £(i,nml)  )  )  /  8.d+00 

ecrest(l)  =  £(i,nm2)  -  y(i,npl) 
y(l,npl)  =  y(l, npl)  +  0.0743801653d+00  *  errest(i) 
go  to  (260,270), jsw 
call  rhs(npl) 
nxt  =  npl 
return 
end 


subroutine  rhs(nxt) 

common  /ham/  t,x( 7, 4 ) , £ (7, 4 ) ,errest ( 7) ,n,h, 101, 102, 103,m, 11, 

+  w,e,sma,u,n4 

double  precision  t,x, £,errest,h, 101, 102, i03,m, ll,w,e,sma,u,n4 
+  Z,a,n4,rl,r,nl,n2,n3,il,i2,l3,b 

dimension  z(7),  r(3),  a(3,3) 


variable  description 


calculate  current  values  o£  inertia 
11  =  101 

12  =  102  +  m  »  (  (0.5*ll*sin(w*t) )**2  ) 

13  =  103  t  m  *  (  (0.5»ll*sin(w*t))**2  ) 

once  againe  convert  variable  names 
z(l)  »  x(l,nxt) 
z(2)  =  x(2,nxt) 


calculate 
a(l, 1) 
a  ( 1, 2 ) 
a  ( 1 ,  3 ) 
a( 2, 1 ) 
a(2,2) 
a  ( 2 ,  3 ) 
a( 3, 1) 
a( 3, 2) 
a  (  3 ,  3 ) 


x( 3,nxt) 
x(4,nxt) 
x(5,nxt) 
x(6,nxt ) 
x(7,nxt) 

;  direction  cosine  matrix  (p5) 


z(l)**2  -  z(  2 )  **2  -  z( 3 ) **2  4-  z(4)**2 
2  *  (  z(l)*z(2)  4  z( 3 ) *z( 4 )  ) 

2  *  (  z(l)*z(3)  -  z( 2 ) *z( 4  )  ) 

2  *  (  z( 1 ) *z( 2 )  -  z( 3 ) *z ( 4 )  ) 

-z(l)**2  4  z(2)**2  -  z( 3 ) **2  4  z( 4 ) *  *  2 

2  *  (  z(2)*z( 3)  +  z(l)*z(4)  ) 

2  *  (  z{ l)*z( 3)  4  z(2)*z(4)  ) 

2  *  (  z(2)*z( 3)  -  z(l)*z(4)  ) 

-z{ 1 ) **2  -  z( 2 ) **2  4  z( 3) **2  4  z(4)**2 


compute  rl  (length  of  orbit)  (pl5) 
rl  =  sma  *  (1  -  e*cos(n4*t)) 

calculate  radius  vector  components  (p8) 
r  (1)  =  a(l,l)*rl 
r ( 2 )  =  a(2,l)*rl 
r ( 3)  =  a(3,l)*rl 

calculate  gravity  gradient  torque  n  (p7) 
b  =  (3*u)  /  ( r 1**5 ) 
nl  =  b  *  (13-12)  *  r ( 2 )  *  r ( 3) 

n2  =  b  *  (11-13)  *  r ( 1 )  *  r ( 3) 

n3  =  b  *  (12-11)  *  r (1)  *  r ( 2) 

write  equations  of  motion  (pp.  6a,  9b) 

f(l,nxt)  =  0.5  *  (  x(2,nxt)*(x(7, nxt ) -a ( 3 , 3 ) *n 4 )  -  x(3,nxt)* 

♦  (x(6,nxt)-a(2,3)*n4)  4  x( 4, nxt ) * ( x( 5, nxt ) -a ( 1, 3 ) *n4 )  ) 

f (2,nxt)  =  0.5  *  (-x(l,nxt)*(x(7,nxt)-a(3,3)*n4)  +  x(3,nxt)» 
t  (x(5,nxt)-a(l,3)*n4)  4-  x( 4, nxt )* (x( 6, nxt ) -a( 2, 3) *n4 )  ) 

f ( 3,nxt )  =  0.5  *  (  x(l,nxt)*(x(6,nxt)-a(2, 3)*n4)  -  x(2,nxt)* 

4  (x(5,nxt)-a(l,3)*n4)  4  x(4,nxt)»(x(7,nxt)-a(3,3)*n4)  ) 

f ( 4 , nxt )  =  0.5  *  (-x(l,nxt)»(x(5,nxt)-a(l,3)»n4)  -  x(2,nxt)* 

4  (x(6,nxt)-a(2,3)*n4)  -  x(3,nxt)*(x(7,nxt)-a(3,3)*n4)  ) 

f (5,nxt)  =  (  nl  +  ( 12-13) *x ( 6,nxt )*x( 7,nxt )  )  /  11 

f ( 6,nxt )  =  (  n2  +  (13-il)*x(7,nxt)*x(5,nxt)  -  ( 0 . 5*w*m*ll**2* 

4  sln(w*t)*cos(v*t) )  *  x(6,nxt)  )  /  12 

f ( 7,nxt )  =  (  n3  +  ( ll-12)»x(5,nxt)*x(6,nxt)  -  ( 0 . 5*v*m*ll**2* 

4  sln(w*t)*cos(w*t) )  *  x(7,nxt)  )  /  13 

return 

end 
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